Parametric

Assessment of Fit of Transitions

We need to derive appropriate functional forms and define respective survival functions. One reason to favor patametric models is that they can be easier to generalize. Several candidate distributions were considered to model transitions, including Weibull (assume a monotonically increasing or decreasing hazard), Log-logistic (non-monotonic hazards), Gompertz (assume a monotonically increasing or decreasing hazard), Log-normal (non-monotonic hazards), Exponential (constant hazard), Gamma, Generalized gamma & Generalized F (more flexible).

The following plots fitted survival curves from each model (colored lines), with the Kaplan-Meier estimate, in red, obtained from an example of Jackson available here, added to the contributions of Wathers & Cutler available here.


#options(warn=-1)

n_iter<-10000

tiempo_antes_fits<-Sys.time()

#Weathers, Brandon and Cutler, Richard Dr., "Comparison of Survival Curves Between Cox Proportional 
#Hazards, Random Forests, and Conditional Inference Forests in Survival Analysis" (2017). All Graduate 
#Plan B and other Reports. 927. 
#https://digitalcommons.usu.edu/gradreports/927 

#<div style="border: 1px solid #ddd; padding: 5px; overflow-y: scroll; height:350px; overflow-x: scroll; width:100%">            
#https://devinincerti.com/2019/01/01/sim-mstate.html

n_trans <- max(trans_matrix, na.rm = TRUE)
fits_wei <- vector(mode = "list", length = n_trans)
fits_weiph <- vector(mode = "list", length = n_trans)
fits_llogis <- vector(mode = "list", length = n_trans)
fits_gomp <- vector(mode = "list", length = n_trans)
fits_logn <- vector(mode = "list", length = n_trans)
fits_exp <- vector(mode = "list", length = n_trans)
fits_gam <- vector(mode = "list", length = n_trans)
fits_ggam <- vector(mode = "list", length = n_trans)
fits_genf <- vector(mode = "list", length = n_trans)
fits_genf_orig <- vector(mode = "list", length = n_trans)
fits_ggam_orig <- vector(mode = "list", length = n_trans)
fits_rp02 <- vector(mode = "list", length = n_trans)
fits_rp03 <- vector(mode = "list", length = n_trans)
fits_rp04 <- vector(mode = "list", length = n_trans)
fits_rp05 <- vector(mode = "list", length = n_trans)
fits_rp06 <- vector(mode = "list", length = n_trans)
fits_rp07 <- vector(mode = "list", length = n_trans)
fits_rp08 <- vector(mode = "list", length = n_trans)
fits_rp09 <- vector(mode = "list", length = n_trans)
fits_rp10 <- vector(mode = "list", length = n_trans)
fits_c_wei <- vector(mode = "list", length = n_trans)
fits_c_weiph <- vector(mode = "list", length = n_trans)
fits_c_llogis <- vector(mode = "list", length = n_trans)
fits_c_gomp <- vector(mode = "list", length = n_trans)
fits_c_logn <- vector(mode = "list", length = n_trans)
fits_c_exp <- vector(mode = "list", length = n_trans)
fits_c_gam <- vector(mode = "list", length = n_trans)
fits_c_ggam <- vector(mode = "list", length = n_trans)
fits_c_genf <- vector(mode = "list", length = n_trans)
fits_c_genf_orig <- vector(mode = "list", length = n_trans)
fits_c_ggam_orig <- vector(mode = "list", length = n_trans)
fits_c_rp02 <- vector(mode = "list", length = n_trans)
fits_c_rp03 <- vector(mode = "list", length = n_trans)
fits_c_rp04 <- vector(mode = "list", length = n_trans)
fits_c_rp05 <- vector(mode = "list", length = n_trans)
fits_c_rp06 <- vector(mode = "list", length = n_trans)
fits_c_rp07 <- vector(mode = "list", length = n_trans)
fits_c_rp08 <- vector(mode = "list", length = n_trans)
fits_c_rp09 <- vector(mode = "list", length = n_trans)
fits_c_rp10 <- vector(mode = "list", length = n_trans)
km.lc <-vector(mode = "list", length = n_trans)

fits_cox<-list()
fits_c_cox<-list()
#"gengamma" Generalized gamma (stable parameterisation)
#"gengamma.orig"    Generalized gamma (original parameterisation)
#"genf" Generalized F (stable parameterisation)
#"genf.orig"    Generalized F (original parameterisation)
#"weibull"  Weibull
#"gamma"    Gamma
#"exp"  Exponential
#"lnorm"    Log-normal
#"gompertz" Gompertz

library(flexsurv)

#Specify the formula
fitform <- Surv(time, status) ~ 1

for (i in 1:n_trans){  
  fits_wei[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibull")
}

for (i in 1:n_trans){  
  fits_weiph[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "weibullph")
}

for (i in 1:n_trans){
  fits_llogis[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "llogis")
}

for (i in 1:n_trans){
  fits_gam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gamma")
}
#LOS ERRORES OCURREN CUANDO NO HAY COVARIABLES
#In (function (q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,  ... :   NaNs produced
#gamma no funcionó: NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced
for (i in 1:n_trans){
  fits_ggam[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma")
}

for (i in 1:n_trans){
  fits_gomp[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "gompertz")
}
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193: valor inicial en 'vmmin' no es finito
for (i in 1:n_trans){
  fits_logn[[i]] <- flexsurvreg(formula=fitform,
                                data = subset(ms_d_match_surv, trans == i),
                                dist = "lnorm")
}

for (i in 1:n_trans){
  fits_exp[[i]] <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "exp")
}

no_attempts <- 20

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "gengamma.orig")
      )
    }
    fits_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
          r <- flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf.orig")
      )
    }
    fits_genf_orig[[i]] <- r
}  
## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

## Warning in flexsurvreg(formula = fitform, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <-  flexsurvreg(formula=fitform,
                               data = subset(ms_d_match_surv, trans == i),
                               dist = "genf")
      )
    }
    fits_genf[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform,k=10,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_rp10[[i]] <- r
}

for (i in 1:n_trans){
km.lc[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
}

transition_label<- plots$title
attr(transition_label,"names")<- plots$trans


get_distinct_hues <- function(ncolor,s=0.5,v=0.95,seed=2125) {
  golden_ratio_conjugate <- 0.618033988749895
  set.seed(seed)
  h <- runif(1)
  H <- vector("numeric",ncolor)
  for(i in seq_len(ncolor)) {
    h <- (h + golden_ratio_conjugate) %% 1
    H[i] <- h
  }
  hsv(H,s=s,v=v)
}
distinct_hues<-get_distinct_hues(20)

layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}
Figure 19. Vissual Assessment of Survival Curves

Figure 19. Vissual Assessment of Survival Curves

endTime <- Sys.time()

paste0("Time in process: ");endTime - tiempo_antes_fits
## [1] "Time in process: "
## Time difference of 3.696909 mins
#23 min aprox.

#For more complicated models, users should specify what covariate values they want summaries for, rather than relying on the default
#</div>
options(warn=0)

if(no_mostrar==1){
jpeg(paste0(dta_path,"_mult_state_ags/exp_surv_int_only.jpg"), height=30, width= 30, res= 600, units = "in")
layout(matrix(1:n_trans, nc = 2, byrow = F))
for (i in 1:n_trans){
plot(km.lc[[i]], col="red", conf.int=F);
  lines(fits_wei[[i]], col=distinct_hues[1], ci=F);
  lines(fits_weiph[[i]], col=distinct_hues[2], ci=F);
  lines(fits_gomp[[i]], col=distinct_hues[3], ci=F);
  lines(fits_llogis[[i]], col=distinct_hues[4], ci=F);#A0A36D
  lines(fits_gam[[i]], col=distinct_hues[5], ci=F);
  lines(fits_ggam[[i]], col=distinct_hues[6], ci=F);
  lines(fits_ggam_orig[[i]], col=distinct_hues[7], ci=F);
  lines(fits_logn[[i]], col=distinct_hues[8], ci=F);
  lines(fits_exp[[i]],col=distinct_hues[9], ci=F);
  lines(fits_genf[[i]],col=distinct_hues[10], ci=F, lty = 2);
  lines(fits_genf_orig[[i]],col=distinct_hues[11], ci=F, lty = 2);  
  lines(fits_rp02[[i]],col=distinct_hues[12], ci=F, lty = 2);
  lines(fits_rp03[[i]],col=distinct_hues[13], ci=F, lty = 2);
  lines(fits_rp04[[i]],col=distinct_hues[14], ci=F, lty = 2);
  lines(fits_rp05[[i]],col=distinct_hues[15], ci=F, lty = 2);
  lines(fits_rp06[[i]],col=distinct_hues[16], ci=F, lty = 2);
  lines(fits_rp07[[i]],col=distinct_hues[17], ci=F, lty = 3);
  lines(fits_rp08[[i]],col=distinct_hues[18], ci=F, lty = 3);
  lines(fits_rp09[[i]],col=distinct_hues[19], ci=F, lty = 3);
  lines(fits_rp10[[i]],col=distinct_hues[20], ci=F, lty = 3)
  
legend("bottomleft", legend = c("Kaplan-Meier","Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (orig)", "Lognormal", "Exponential", "Generalized F","Generalized F (orig)","RP2","RP3","RP4","RP5", "RP6","RP7", "RP8", "RP9", "RP10"), col = 
         c("red",distinct_hues), 
       title = "Distributions", cex = .95, bty = "n", lty=c(rep(1,10),rep(2,6),rep(3,4)),lwd=3)# lty = 1:2, 
title(main=transition_label[[i]])
}
  dev.off()
}

The following transitions deviated considerably from the models chosen. Consider specify transition parameters manually or check whether to transform the measure of time into an exponential or logarithmic time to event:

  • Readmission2 to DWCA3
  • Readmission2 to TD3
  • Readmission to DWCA2
  • Readmission to TD2
  • Admission to Readmission
  • Admission to DWCA
  • Admission to TD


#
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
fit_flexsurvreg0<-data.frame()
fitted_flexsurvreg0<-data.frame()

dists_no_covs_10s<-cbind.data.frame(covs=c(rep("fits_",(20)*n_trans)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans),
              dist=c("weibull", "weibullph", "llogis", "gamma", "gengamma","gengamma.orig", "gompertz", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans),
              trans=rep(1:n_trans, each=20))

for (i in 1:nrow(dists_no_covs_10s)){  #
  
  cat(paste0("#### Flexible Survival Model (w/ covs): ",
               dists_no_covs_10s[i,"formal"], "; transition: ",dists_no_covs_10s[i,"trans"], "\n \n"))  
    
  model<-paste0("fits_",dists_no_covs_10s[i,"model"])
  
  if(!is.null(get(model)[[dists_no_covs_10s[i,"trans"]]])){  
     fitted_flexsurvreg0<- rbind(fitted_flexsurvreg0,cbind.data.frame(dist=rep(dists_no_covs_10s[i,"formal"],), 
                            trans=rep(dists_no_covs_10s[i,"trans"],),
                            data.table::data.table(summary(get(model)[[dists_no_covs_10s[i,"trans"]]], 
                                                           #newdata= newdat2a, 
                                                           #newtime=unique(fitted_km$time), 
                                                           type = "survival", 
                                                           tidy=T)))) 
    # Generate fit indices
    fit_flexsurvreg0<-rbind(fit_flexsurvreg0,
       cbind(dist= dists_no_covs_10s[i,"formal"],
             transition=dists_no_covs_10s[i,"trans"],
             fitstats.flexsurvreg(get(model)[[dists_no_covs_10s[i,"trans"]]])))
    #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
    } else {
    cat(paste0("The model that assumed a ",dists_no_covs_10s[i,"formal"]," distribution for the transition number ",dists_no_covs_10s[i,"trans"]," did not converge \n\n"))
    }
}
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## The model that assumed a Gompertz distribution for the transition number 2 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## The model that assumed a Gompertz distribution for the transition number 3 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## The model that assumed a Gompertz distribution for the transition number 4 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg0 %>% 
    dplyr::group_by(transition) %>% 
    summarise(mean(ucl,na.rm=T))
  }
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

km.lc0<-list()
fitted_km0<-data.frame()

for (i in 1:n_trans){
  km.lc0[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i))
  fitted_km0<-rbind(fitted_km0,cbind.data.frame(trans=i,fortify(km.lc0[[i]])))
}

#Calculate error
#c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma","Generalized gamma", "Lognormal", "Exponential", "Generalized F")

fitted_flexsurvreg_binned_mix0<-
data.frame(fitted_flexsurvreg0)[,c("dist","trans","time","est","lcl","ucl")] %>% 
  dplyr::left_join(fitted_km0, by=c("trans","time")) %>% 
#there are many observations that was not available due to empirical missing data
#dplyr::filter(!is.na(surv))
  dplyr::group_by(dist,trans) %>% 
  #dplyr::mutate(est= ifelse(is.na(est), mean(est, na.rm=TRUE), est)) %>% 
  #dplyr::mutate(surv= ifelse(is.na(surv), mean(surv, na.rm=TRUE), surv)) %>% 
  tidyr::fill(surv,.direction="updown") %>% 
  tidyr::fill(est,.direction="updown") %>% 
  ungroup()

db_for_apply_rmse0<-
  cbind.data.frame(
      dist=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F (orignial)", paste0("RP0",2:9),"RP10"),n_trans),
      trans=rep(c(1:n_trans),each=20*2))
   
rmse_comp_fits0<- data.frame()
for(i in 1:nrow(db_for_apply_rmse0)){
  rmse<- Metrics::rmse(subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$est,
              subset(fitted_flexsurvreg_binned_mix0, dist==db_for_apply_rmse0[i,"dist"] & 
                       trans==db_for_apply_rmse0[i,"trans"])$surv)

  rmse_comp_fits0<- rbind(rmse_comp_fits0,cbind(dist=db_for_apply_rmse0[i,"dist"],
                                                  residential=db_for_apply_rmse0[i,"residential"],
                                                  trans=db_for_apply_rmse0[i,"trans"],
                                                  rmse=rmse))
}

rmse_comp_fits0_filter<-
  rmse_comp_fits0 %>% 
    dplyr::filter(rmse!="NaN") %>%  
      dplyr::arrange(trans,dist, rmse)%>%
      dplyr::mutate(rmse=round(as.numeric(rmse),4))

setting_type_label<- c(`0`="Outpatient",`1`="Residential") 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_


We checked the AICs, cumulative relative error (differences between the Kaplan-Meier curve, and the observed behavior of each distribution).


#load("C:/Users/CISS Fondecyt/OneDrive/Escritorio/mult_state_ags.RData")

c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_binned_mix0_plot<-
fitted_flexsurvreg_binned_mix0 %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% group_by(dist,trans) %>% dplyr::mutate(cumsum_error=cumsum(`abs((est - surv)/surv)`)) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist),size=1)+
    geom_line(aes(time, surv), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
                         #c("black",brewer.pal(n = 9, name = 'Paired')))+
                         #c("#112A60","#085754","#D3A347","#4F3C91","red","#112A60","#085754","#8F630D","#251363")) +
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    sjPlot::theme_sjplot2()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
  scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  #ylim(0,1.25)+
  #theme(axis.text.x = element_blank(), 
  #      panel.grid.major = element_blank(), 
  #      panel.grid.minor = element_blank()) +
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)
p1<-
fitted_flexsurvreg_binned_mix0_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.03, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))

p1

Figure 20. Vissual Assessment of Survival Curves (w/o covars)

  #htmlwidgets::saveWidget(as_widget(partial_bundle(p1, local=T)), "test1.html")

#https://github.com/ropensci/plotly/issues/1586

                 #%>% 
  #layout(annotations = list(
  #              list(x = 0.5 , y = -0.03, text = "Survival",
  #                   font = list(color = "black",size = 15),
  #                   textangle = 0,
  #                   showarrow = F, xref='paper', yref='paper', size=48,
  #                   automargin = T)
   #             ))
n_trans2 <- max(trans_matrix, na.rm = T)

fit_flexurvreg0_kable<-
fit_flexsurvreg0 %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb, dist, Df, n2ll, AIC, AICc, BIC) %>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))
  
 
fit_flexurvreg0_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 11. Fit indices of the adjusted survival analyses (intercept-only)"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 
Table 11. Fit indices of the adjusted survival analyses (intercept-only)
Transition Distribution DF Negative 2 Log Likelihood AIC AICc BIC
01)Admission to Readmission Generalized F(original) 4 NA NA NA NA
01)Admission to Readmission RP06 8 121036.84 121052.84 121052.85 121117.16
01)Admission to Readmission RP04 6 121041.85 121053.85 121053.86 121102.09
01)Admission to Readmission RP08 10 121033.99 121053.99 121054.00 121134.39
01)Admission to Readmission Generalized F 4 121046.51 121054.51 121054.51 121086.67
01)Admission to Readmission RP10 12 121031.52 121055.52 121055.53 121151.99
01)Admission to Readmission RP09 11 121033.96 121055.96 121055.97 121144.39
01)Admission to Readmission RP07 9 121038.72 121056.72 121056.73 121129.08
01)Admission to Readmission RP05 7 121043.64 121057.64 121057.65 121113.92
01)Admission to Readmission RP03 5 121058.64 121068.64 121068.65 121108.84
01)Admission to Readmission RP02 4 121068.72 121076.72 121076.72 121108.88
01)Admission to Readmission Generalized gamma 3 121291.72 121297.72 121297.72 121321.84
01)Admission to Readmission Lognormal 2 121619.49 121623.49 121623.49 121639.57
01)Admission to Readmission Generalized gamma (original) 3 121759.04 121765.04 121765.05 121789.16
01)Admission to Readmission Log-logistic 2 122122.85 122126.85 122126.85 122142.93
01)Admission to Readmission Weibull (AFT) 2 122464.79 122468.79 122468.79 122484.87
01)Admission to Readmission Weibull (PH) 2 122464.79 122468.79 122468.79 122484.87
01)Admission to Readmission Gamma 2 122578.49 122582.49 122582.49 122598.57
01)Admission to Readmission Exponential 1 122964.89 122966.89 122966.89 122974.93
02)Readmission to Readmission2 Generalized F 4 36551.44 36559.44 36559.45 36586.50
02)Readmission to Readmission2 Generalized F(original) 4 36551.44 36559.44 36559.45 36586.50
02)Readmission to Readmission2 RP03 5 36551.19 36561.19 36561.20 36595.01
02)Readmission to Readmission2 RP04 6 36549.46 36561.46 36561.48 36602.05
02)Readmission to Readmission2 RP05 7 36547.95 36561.95 36561.96 36609.29
02)Readmission to Readmission2 RP08 10 36543.34 36563.34 36563.37 36630.97
02)Readmission to Readmission2 RP09 11 36541.54 36563.54 36563.58 36637.94
02)Readmission to Readmission2 RP06 8 36547.72 36563.72 36563.74 36617.83
02)Readmission to Readmission2 RP07 9 36545.91 36563.91 36563.94 36624.78
02)Readmission to Readmission2 RP10 12 36541.57 36565.57 36565.62 36646.73
02)Readmission to Readmission2 RP02 4 36561.33 36569.33 36569.34 36596.39
02)Readmission to Readmission2 Generalized gamma 3 36651.53 36657.53 36657.54 36677.82
02)Readmission to Readmission2 Lognormal 2 36676.37 36680.37 36680.37 36693.90
02)Readmission to Readmission2 Generalized gamma (original) 3 36704.13 36710.13 36710.13 36730.42
02)Readmission to Readmission2 Log-logistic 2 36786.43 36790.43 36790.43 36803.96
02)Readmission to Readmission2 Weibull (AFT) 2 36905.75 36909.75 36909.76 36923.28
02)Readmission to Readmission2 Weibull (PH) 2 36905.75 36909.75 36909.76 36923.28
02)Readmission to Readmission2 Gamma 2 36931.00 36935.00 36935.00 36948.53
02)Readmission to Readmission2 Exponential 1 36978.75 36980.75 36980.75 36987.51
03)Readmission2 to Readmission3 Generalized F(original) 4 NA NA NA NA
03)Readmission2 to Readmission3 Generalized F 4 11510.25 11518.25 11518.27 11540.68
03)Readmission2 to Readmission3 RP02 4 11511.22 11519.22 11519.24 11541.64
03)Readmission2 to Readmission3 RP03 5 11511.36 11521.36 11521.39 11549.39
03)Readmission2 to Readmission3 RP04 6 11510.57 11522.57 11522.61 11556.21
03)Readmission2 to Readmission3 RP05 7 11510.45 11524.45 11524.51 11563.70
03)Readmission2 to Readmission3 RP06 8 11509.57 11525.57 11525.64 11570.42
03)Readmission2 to Readmission3 RP07 9 11509.11 11527.11 11527.20 11577.56
03)Readmission2 to Readmission3 RP08 10 11507.51 11527.51 11527.62 11583.57
03)Readmission2 to Readmission3 RP09 11 11506.53 11528.53 11528.66 11590.20
03)Readmission2 to Readmission3 RP10 12 11505.30 11529.30 11529.46 11596.58
03)Readmission2 to Readmission3 Generalized gamma 3 11532.03 11538.03 11538.04 11554.85
03)Readmission2 to Readmission3 Lognormal 2 11541.52 11545.52 11545.53 11556.73
03)Readmission2 to Readmission3 Generalized gamma (original) 3 11552.33 11558.33 11558.34 11575.15
03)Readmission2 to Readmission3 Log-logistic 2 11578.52 11582.52 11582.53 11593.73
03)Readmission2 to Readmission3 Weibull (AFT) 2 11621.91 11625.91 11625.92 11637.12
03)Readmission2 to Readmission3 Weibull (PH) 2 11621.91 11625.91 11625.92 11637.12
03)Readmission2 to Readmission3 Gamma 2 11626.91 11630.91 11630.91 11642.12
03)Readmission2 to Readmission3 Exponential 1 11630.65 11632.65 11632.66 11638.26
04)Readmission3 to Readmission4 Generalized F(original) 4 3607.60 3615.60 3615.67 3633.49
04)Readmission3 to Readmission4 Generalized F 4 3607.67 3615.67 3615.73 3633.56
04)Readmission3 to Readmission4 RP02 4 3608.19 3616.19 3616.25 3634.08
04)Readmission3 to Readmission4 RP03 5 3608.33 3618.33 3618.42 3640.69
04)Readmission3 to Readmission4 RP04 6 3607.92 3619.92 3620.05 3646.75
04)Readmission3 to Readmission4 RP05 7 3606.91 3620.91 3621.08 3652.21
04)Readmission3 to Readmission4 RP06 8 3606.40 3622.40 3622.63 3658.18
04)Readmission3 to Readmission4 Lognormal 2 3619.93 3623.93 3623.95 3632.88
04)Readmission3 to Readmission4 RP07 9 3606.32 3624.32 3624.60 3664.57
04)Readmission3 to Readmission4 Generalized gamma 3 3619.45 3625.45 3625.48 3638.86
04)Readmission3 to Readmission4 RP08 10 3605.47 3625.47 3625.82 3670.20
04)Readmission3 to Readmission4 RP09 11 3605.77 3627.77 3628.18 3676.96
04)Readmission3 to Readmission4 Generalized gamma (original) 3 3622.23 3628.23 3628.27 3641.65
04)Readmission3 to Readmission4 RP10 12 3605.95 3629.95 3630.44 3683.61
04)Readmission3 to Readmission4 Log-logistic 2 3627.00 3631.00 3631.02 3639.95
04)Readmission3 to Readmission4 Exponential 1 3641.53 3643.53 3643.54 3648.00
04)Readmission3 to Readmission4 Weibull (AFT) 2 3640.15 3644.15 3644.17 3653.10
04)Readmission3 to Readmission4 Weibull (PH) 2 3640.15 3644.15 3644.17 3653.10
04)Readmission3 to Readmission4 Gamma 2 3641.13 3645.13 3645.14 3654.07


rmse_plot0<-
rmse_comp_fits0_filter %>% 
  dplyr::mutate(trans=as.numeric(trans)) %>% 
  dplyr::mutate(trans2=as.numeric(trans)) %>% 
    dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("trans"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",trans),")",transition_label)) %>% 
        
   ggplot()+
  geom_bar(aes(x=dist, y=rmse), position="dodge", stat="identity", alpha=0.4)+
  sjPlot::theme_sjplot2()+
  #geom_errorbar(aes(x=t, ymin=L, ymax=U, color=program), position="dodge", stat="identity", width=0.4,  alpha=0.9, size=1.3)+
  facet_wrap(.~trans_w_nmb, ncol=3, scales="free_y", dir="h") + 
  xlab("") + 
  ylab("")+
  #ylab("State occupancy probabilities") + 
  theme_minimal()+
  theme(legend.position="bottom")+
    theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        panel.border = element_blank(),
          legend.box.background = element_blank())+
  theme(axis.text.x = element_text(angle = 90,hjust=0.95,vjust=0.2))#, vjust=0, hjust=0.3

rmse_plot0
Figure 21. Suvival RMSEs, Ten-states Model (w/o covars)

Figure 21. Suvival RMSEs, Ten-states Model (w/o covars)

#install.packages("survHE")
library(survHE)
#First, defines the vector of models to be used
mods <- c("weibull", "weibullph", "llogis", "gamma", "gengamma", "gompertz", "lnorm" ,"exp")#, "genf")

# And then runs the models using MLE via flexsurv
#m2 <- fit.models(formula = formula, data = data, distr = c("exp","wei","lno","llo"), method="inla")
#m3 <- fit.models(formula = formula, data = subset(ms_d_match_surv, trans == 1), distr = mods, method="hmc")
#Error in model.matrix(formula, data)[(mf %>% filter(event == 1))$ID, ] : 
#  subíndice fuera de  los límites

list_aics<-list()
for (i in 1:n_trans){
  m1 <- fit.models(formula = fitform, data = subset(ms_d_match_surv, trans == i), distr = mods)
  model.fit.plot(MLE=m1, stacked = T, #HMC=m3, 
               models = c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Lognormal", "Exponential", "Generalized F"),
               #lab.profile = c("Outpatient","Residential")
               name_legend = "Inferential method")+ 
  ggtitle("Model comparison based on AIC")+
  labs(subtitle=plots[i,"title"])
}


tiempo_antes_fits2<-Sys.time()

newtime0 = seq(from=round(min(ms_d_match_surv$time),0), to= round(max(ms_d_match_surv$time),0))

#_#_#_#_#_#_#_#_#_#_#_#_#_
#covariates

#Specify the formula
fitform2 <- Surv(time, status) ~  factor(tipo_de_plan_res_1)


kernel_haz_est2a<-list()
kernel_haz_est2b<-list()
kernel_haz2a<-list()
kernel_haz2b<-list()
for (i in 1:n_trans){
library("muhaz")
kernel_haz_est2a[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==1),"status"])
kernel_haz2a[[i]] <- data.table(time = kernel_haz_est2a[[i]]$est.grid,
                         est = kernel_haz_est2a[[i]]$haz.est,
                         dist = "Kernel density")
kernel_haz_est2b[[i]] <- muhaz(ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"time"],
                        ms_d_match_surv[which(ms_d_match_surv$trans==i &
                        ms_d_match_surv$tipo_de_plan_res_1==0),"status"])
kernel_haz2b[[i]] <- data.table(time = kernel_haz_est2b[[i]]$est.grid,
                         est = kernel_haz_est2b[[i]]$haz.est,
                         dist = "Kernel density")
}

haz_int_only2<-
  rbind(cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2a[[i]])),
              tipo_de_plan_res_1=rep(1,nrow(kernel_haz2a[[i]])),
      rbindlist(kernel_haz2a)),
      cbind(trans=rep(1:n_trans,each=nrow(kernel_haz2b[[i]])),
            tipo_de_plan_res_1=rep(0,nrow(kernel_haz2b[[i]])),
      rbindlist(kernel_haz2b)))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

n_trans2 <- max(trans_matrix, na.rm = T)

dists_w_covs_10s<-cbind.data.frame(covs=c(rep("fits_c_",(20)*n_trans2)),
              formal=rep(c("Weibull (AFT)", "Weibull (PH)", "Gompertz", "Log-logistic", "Gamma",
                "Generalized gamma", "Generalized gamma (original)", "Lognormal", "Exponential", "Generalized F", "Generalized F(original)", paste0("RP0",2:9),"RP10"),1*n_trans2),
              dist=c("weibull", "weibullph", "llogis", "gamma", "gengamma","gengamma.orig", "gompertz", "lnorm", "exp", "genf","genf.orig",rep("no dist",9)),
              model=rep(c("wei", "weiph", "gomp", "llogis", "gam","ggam", "ggam_orig", "logn", "exp", "genf","genf_orig", paste0("rp0",2:9),"rp10"),1*n_trans2),
              trans=rep(1:n_trans2, each=20))

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
no_attempts <- 30

km.lc2a<-list()
km.lc2b<-list()
fits_cox2<-list()
fits_c_cox2<-list()

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibull")
    )
  } 
  fits_c_wei[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "weibullph")
    )
  } 
  fits_c_weiph[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "llogis")
    )
  } 
  fits_c_llogis[[i]] <- r
}

for (i in 1:n_trans2){
  r <- NULL
  attempt <- 0
  while( is.null(r) && attempt <= no_attempts ) {
    attempt <- attempt + 1
    try(
      r <- flexsurvreg(formula=fitform2,
                                 data = subset(ms_d_match_surv, trans == i),
                                 dist = "gamma")
    )
  } 
  fits_c_gam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "gengamma")
      )
    }
    fits_c_ggam[[i]] <- r
}

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i), 
                                   dist = "gompertz")
      )
    }
    fits_c_gomp[[i]] <- r
}  
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
## Error in optim(method = "BFGS", par = c(shape = 0.001, rate = -7.33378535554193,  : 
##   valor inicial en 'vmmin' no es finito
for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "lnorm")
      )
    }
    fits_c_logn[[i]] <- r
}  

for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "exp")
      )
    }
    fits_c_exp[[i]] <- r
}  


for (i in 1:n_trans2){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvreg(formula=fitform2,
                                   data = subset(ms_d_match_surv, trans == i),
                                   dist = "genf")
      )
    }
    fits_c_genf[[i]] <- r
}  


for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "gengamma.orig")
        )
    }
    fits_c_ggam_orig[[i]] <- r
}  

for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
        attempt <- attempt + 1
        try(
            r <- flexsurvreg(formula=fitform2,
                             data = subset(ms_d_match_surv, trans == i),
                             dist = "genf.orig")
        )
    }
    fits_c_genf_orig[[i]] <- r
}
## Warning in flexsurvreg(formula = fitform2, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.

## Warning in flexsurvreg(formula = fitform2, data = subset(ms_d_match_surv, :
## Optimisation has probably not converged to the maximum likelihood - Hessian is
## not positive definite.
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=2,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp02[[i]] <- r
}  
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=3,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp03[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=4,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp04[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=5,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp05[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=6,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp06[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=7,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp07[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=8,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp08[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=9,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp09[[i]] <- r
}
for (i in 1:n_trans){
    r <- NULL
    attempt <- 0
    while( is.null(r) && attempt <= no_attempts ) {
      attempt <- attempt + 1
      try(
        r <- flexsurvspline(formula=fitform2,k=10,
                               data = subset(ms_d_match_surv, trans == i))
      )
    }
    fits_c_rp10[[i]] <- r
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

fitted_km<-data.frame()
for (i in 1:n_trans2){
km.lc2a[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==1))
km.lc2b[[i]] <- survfit(formula= fitform, data = subset(ms_d_match_surv, trans == i & tipo_de_plan_res_1==0))

fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(1,), fortify(km.lc2a[[i]])))
fitted_km<-rbind(fitted_km,cbind.data.frame(trans=i,residential=rep(0,), fortify(km.lc2b[[i]])))
}

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

newdat2a <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(1,1))))
newdat2b <- data.table::data.table(tipo_de_plan_res_1= factor(c(rep(0,1))))

fitted_flexsurvreg<-data.frame()
fit_flexsurvreg<-data.frame()

for (i in 1:nrow(dists_w_covs_10s)){  #
  
cat(paste0("#### Flexible Survival Model (w/ covs): ",
             dists_w_covs_10s[i,"formal"], "; transition: ",dists_w_covs_10s[i,"trans"], "\n \n"))  
  
model<-paste0("fits_c_",dists_w_covs_10s[i,"model"])

if(!is.null(get(model)[[dists_w_covs_10s[i,"trans"]]])){  
  #Generate databases
 fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(1,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2a, newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
  fitted_flexsurvreg<- rbind(fitted_flexsurvreg,cbind.data.frame(dist=rep(dists_w_covs_10s[i,"formal"],), 
                            trans=rep(dists_w_covs_10s[i,"trans"],),
                            residential=rep(0,),
                            data.table::data.table(summary(get(model)[[dists_w_covs_10s[i,"trans"]]], newdata= newdat2b,  newtime=unique(fitted_km$time), type = "survival", tidy=T)))) 
 #t=newtime0, 
 
  # Generate fit indices
  fit_flexsurvreg<-rbind(fit_flexsurvreg,
     cbind(dist= dists_w_covs_10s[i,"formal"],
           transition=dists_w_covs_10s[i,"trans"],
           fitstats.flexsurvreg(get(model)[[dists_w_covs_10s[i,"trans"]]])))
  #the BIC may not be appropriate if none of the candidate models are considered to be close to the ‘true’ model.     
  } else {
  cat(paste0("The model that assumed a ",dists_w_covs_10s[i,"formal"]," distribution for the transition number ",dists_w_covs_10s[i,"trans"]," did not converge \n\n"))
  }
}
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 1
##  
## The model that assumed a Gompertz distribution for the transition number 1 did not converge 
## 
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 1
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 1
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 2
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 2
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 3
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 3
##  
## #### Flexible Survival Model (w/ covs): Weibull (AFT); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Weibull (PH); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gompertz; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Log-logistic; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized gamma (original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): Lognormal; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Exponential; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F; transition: 4
##  
## #### Flexible Survival Model (w/ covs): Generalized F(original); transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP02; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP03; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP04; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP05; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP06; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP07; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP08; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP09; transition: 4
##  
## #### Flexible Survival Model (w/ covs): RP10; transition: 4
## 
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(no_mostrar==1){
  fit_flexsurvreg %>% 
    dplyr::group_by(trans) %>% 
    summarise(mean(ucl,na.rm=T))
  }

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_


#load("C:/Users/CISS Fondecyt/OneDrive/Escritorio/mult_state_ags.RData")

c26 <- c(
  "dodgerblue2", # red
  "green4",
  "#6A3D9A", # purple
  "#FF7F00", # orange
  "gray16", "gold1",
  "skyblue2", "#FB9A99", # lt pink
  "palegreen2",
  "#CAB2D6", # lt purple
  "#FDBF6F", # lt orange
  "gray70", "khaki2",
  "maroon", "orchid1", "deeppink1", "blue1", "steelblue4",
  "darkturquoise", "green1", "yellow4", "yellow3",
  "darkorange4", "brown", "gray40")

fitted_flexsurvreg_plot<-
fitted_flexsurvreg %>%  
  dplyr::left_join(fitted_km, by=c("time", "trans", "residential")) %>% 
  group_by(trans,dist,residential) %>% 
  tidyr::fill(surv,.direction="down") %>% 
  ungroup() %>% 
  dplyr::mutate(abs((est-surv)/surv)) %>% 
  group_by(trans,dist,residential) %>% 
  dplyr::mutate(cumsum_error=cumsum(replace_na(`abs((est - surv)/surv)`, 0))) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(text=paste0("Distribution: ",dist,"<br>survival: ",sprintf("%1.2f",est),"<br>Time (in years): ",round(time/365.25,2), "<br>CIs: ",ifelse(!is.na(lcl),"CIs","no CIs"),"<br>Cumsum(abs((est-surv)/surv))= ",round(cumsum_error,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>% 
  dplyr::mutate(text2=paste0("Distribution: KM","<br>survival: ",sprintf("%1.2f",surv),"<br>Time (in years): ",round(time/365.25,2),"<br>Residential: ",ifelse(residential==1,"Yes","No"))) %>%
ggplot()+
    geom_line(aes(time, est, color=dist, linetype=factor(residential)),size=1)+
    geom_line(aes(time, surv, linetype=factor(residential)), color="red",size=1)+
    scale_color_manual(name="Distributions", values = c26[1:20])+
    facet_wrap(.~trans,labeller = labeller(trans = transition_label), ncol=3, scales = "free_y")+
    sjPlot::theme_sjplot2()+
    theme(legend.position="bottom",
          strip.background = element_rect(fill = "white", colour = "white"))+
 # scale_x_continuous(breaks = seq(0, 365.25*12, 2*365.25), labels=seq(0,12,2))+
  labs(y="",x="")+
  theme(legend.position='none')

m <- list(
  l = 80*1.05,
  r = 80,
  b = 80,
  t = 80*1.05,
  pad = 4
)

p2<-
fitted_flexsurvreg_plot$data %>% 
  dplyr::mutate(years=round(time/365.25,0)) %>% 
  dplyr::mutate(trans2=trans) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans2=1:n_trans2), by="trans2") %>% 
group_by(trans) %>%
group_map(~ plot_ly(data=., x = ~time, y = ~est, color = ~dist, type = "scatter", mode="lines", linetype = ~residential, text=~text) %>% 
            add_trace(y = ~surv, type = 'scatter', mode = 'lines', line = list(color = 'red', width = 1, linetype = ~residential), text=~text2)%>% 
  layout(annotations = list(list(x = 0.5 , y = 1.03, text = ~unique(transition_label), showarrow = F, xref='paper', yref='paper'))) %>% 
  layout(
    xaxis = list(
      ticktext = list(1, 3, 5, 7, 9, 11), 
      textangle = 0,
      tickvals = list(365.25*1, 365.25*3, 365.25*5, 365.25*7, 365.25*9, 365.25*11),
      tickmode = "array"
  )), keep=TRUE) %>%    
subplot(nrows = 3, shareX = T, shareY=T, titleX = F, titleY = F, margin = .05)%>% layout(showlegend = F) %>% #, margin = 0.05) %>% 
  layout(annotations = list(
                list(x = -0.1, y = 0.5, text = "Survival",
                     font = list(color = "black",size = 15),
                     textangle = 270,
                     showarrow = F, xref='paper', yref='paper', size=48,margin =m)))
## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

## Warning in RColorBrewer::brewer.pal(N, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors
#p2

  htmlwidgets::saveWidget(as_widget(p2), "test2_jun.html")
#https://docs.ropensci.org/plotly/reference/partial_bundle.html

File is available in this link.

fit_flexurvreg_kable<-
fit_flexsurvreg %>% 
  dplyr::arrange(transition, AIC) %>% 
  dplyr::left_join(cbind.data.frame(transition_label,trans_nmb=1:n_trans2),by=c("transition"="trans_nmb")) %>% 
  dplyr::mutate(trans_w_nmb=paste0(sprintf("%02d",transition),")",transition_label)) %>% 
  dplyr::select(trans_w_nmb,dist,Df,n2ll,AIC,AICc,BIC)%>%
  dplyr::mutate_at(vars(n2ll, AIC, AICc, BIC),~ifelse(abs(.)>1e6,NA,format(round(.,2), scientific = FALSE)))  
  
fit_flexurvreg_kable %>%   
      knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption = paste0("Table 12. Fit indices of the adjusted survival analyses"),
               col.names = c("Transition","Distribution", "DF","Negative 2 Log Likelihood","AIC","AICc","BIC"),
               align =c("l",rep('c', 101))) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12) %>%
  #kableExtra::add_footnote("Note. NA= Null values", notation="none") %>% 
  kableExtra::scroll_box(width = "100%", height = "375px") 
Table 12. Fit indices of the adjusted survival analyses
Transition Distribution DF Negative 2 Log Likelihood AIC AICc BIC
01)Admission to Readmission Generalized F(original) 5 NA NA NA NA
01)Admission to Readmission Generalized F 5 120769.03 120779.03 120779.03 120819.23
01)Admission to Readmission RP06 9 120779.16 120797.16 120797.17 120869.52
01)Admission to Readmission RP04 7 120784.24 120798.24 120798.25 120854.52
01)Admission to Readmission RP08 11 120776.37 120798.37 120798.38 120886.81
01)Admission to Readmission RP10 13 120773.93 120799.93 120799.95 120904.45
01)Admission to Readmission RP09 12 120776.36 120800.36 120800.37 120896.83
01)Admission to Readmission RP07 10 120781.05 120801.05 120801.06 120881.44
01)Admission to Readmission RP05 8 120786.03 120802.03 120802.04 120866.35
01)Admission to Readmission RP03 6 120800.94 120812.94 120812.95 120861.18
01)Admission to Readmission RP02 5 120811.14 120821.14 120821.15 120861.34
01)Admission to Readmission Generalized gamma 4 120965.59 120973.59 120973.59 121005.74
01)Admission to Readmission Lognormal 3 121321.19 121327.19 121327.19 121351.31
01)Admission to Readmission Generalized gamma (original) 4 121472.61 121480.61 121480.61 121512.76
01)Admission to Readmission Log-logistic 3 121840.98 121846.98 121846.98 121871.10
01)Admission to Readmission Weibull (AFT) 3 122215.75 122221.75 122221.76 122245.87
01)Admission to Readmission Weibull (PH) 3 122215.75 122221.75 122221.76 122245.87
01)Admission to Readmission Gamma 3 122337.40 122343.40 122343.40 122367.52
01)Admission to Readmission Exponential 2 122716.58 122720.58 122720.59 122736.66
02)Readmission to Readmission2 Generalized F 5 36551.26 36561.26 36561.27 36595.07
02)Readmission to Readmission2 Generalized F(original) 5 36551.26 36561.26 36561.27 36595.07
02)Readmission to Readmission2 RP03 6 36550.53 36562.53 36562.55 36603.12
02)Readmission to Readmission2 RP04 7 36548.81 36562.81 36562.83 36610.15
02)Readmission to Readmission2 RP05 8 36547.29 36563.29 36563.31 36617.40
02)Readmission to Readmission2 RP08 11 36542.68 36564.68 36564.72 36639.08
02)Readmission to Readmission2 RP09 12 36540.81 36564.81 36564.85 36645.97
02)Readmission to Readmission2 RP06 9 36547.06 36565.06 36565.09 36625.93
02)Readmission to Readmission2 RP07 10 36545.25 36565.25 36565.29 36632.89
02)Readmission to Readmission2 RP10 13 36540.91 36566.91 36566.97 36654.84
02)Readmission to Readmission2 RP02 5 36560.67 36570.67 36570.68 36604.48
02)Readmission to Readmission2 Generalized gamma 4 36650.99 36658.99 36658.99 36686.04
02)Readmission to Readmission2 Lognormal 3 36675.82 36681.82 36681.82 36702.11
02)Readmission to Readmission2 Gompertz 3 36697.35 36703.35 36703.35 36723.64
02)Readmission to Readmission2 Generalized gamma (original) 4 36707.55 36715.55 36715.55 36742.60
02)Readmission to Readmission2 Log-logistic 3 36785.86 36791.86 36791.86 36812.15
02)Readmission to Readmission2 Weibull (AFT) 3 36905.33 36911.33 36911.33 36931.62
02)Readmission to Readmission2 Weibull (PH) 3 36905.33 36911.33 36911.33 36931.62
02)Readmission to Readmission2 Gamma 3 36930.62 36936.62 36936.62 36956.91
02)Readmission to Readmission2 Exponential 2 36978.44 36982.44 36982.44 36995.97
03)Readmission2 to Readmission3 Generalized F(original) 5 11508.11 11518.11 11518.14 11546.14
03)Readmission2 to Readmission3 Generalized F 5 11508.11 11518.11 11518.14 11546.14
03)Readmission2 to Readmission3 RP02 5 11511.21 11521.21 11521.24 11549.24
03)Readmission2 to Readmission3 RP03 6 11511.36 11523.36 11523.40 11557.00
03)Readmission2 to Readmission3 RP04 7 11510.56 11524.56 11524.62 11563.81
03)Readmission2 to Readmission3 RP05 8 11510.45 11526.45 11526.52 11571.30
03)Readmission2 to Readmission3 RP06 9 11509.57 11527.57 11527.66 11578.03
03)Readmission2 to Readmission3 RP07 10 11509.10 11529.10 11529.21 11585.17
03)Readmission2 to Readmission3 RP08 11 11507.51 11529.51 11529.64 11591.18
03)Readmission2 to Readmission3 RP09 12 11506.53 11530.53 11530.68 11597.80
03)Readmission2 to Readmission3 RP10 13 11505.29 11531.29 11531.47 11604.18
03)Readmission2 to Readmission3 Generalized gamma 4 11530.60 11538.60 11538.62 11561.02
03)Readmission2 to Readmission3 Lognormal 3 11540.98 11546.98 11546.99 11563.80
03)Readmission2 to Readmission3 Generalized gamma (original) 4 11552.49 11560.49 11560.51 11582.92
03)Readmission2 to Readmission3 Gompertz 3 11568.15 11574.15 11574.16 11590.97
03)Readmission2 to Readmission3 Log-logistic 3 11578.40 11584.40 11584.42 11601.22
03)Readmission2 to Readmission3 Weibull (AFT) 3 11621.90 11627.90 11627.91 11644.72
03)Readmission2 to Readmission3 Weibull (PH) 3 11621.90 11627.90 11627.91 11644.72
03)Readmission2 to Readmission3 Gamma 3 11626.88 11632.88 11632.90 11649.70
03)Readmission2 to Readmission3 Exponential 2 11630.63 11634.63 11634.63 11645.84
04)Readmission3 to Readmission4 Generalized F(original) 5 3607.24 3617.24 3617.34 3639.61
04)Readmission3 to Readmission4 RP02 5 3607.36 3617.36 3617.45 3639.72
04)Readmission3 to Readmission4 Generalized F 5 3607.61 3617.61 3617.70 3639.97
04)Readmission3 to Readmission4 RP03 6 3607.49 3619.49 3619.62 3646.32
04)Readmission3 to Readmission4 RP04 7 3607.08 3621.08 3621.26 3652.39
04)Readmission3 to Readmission4 RP05 8 3606.09 3622.09 3622.32 3657.87
04)Readmission3 to Readmission4 RP06 9 3605.59 3623.59 3623.88 3663.84
04)Readmission3 to Readmission4 RP07 10 3605.51 3625.51 3625.86 3670.24
04)Readmission3 to Readmission4 Lognormal 3 3619.70 3625.70 3625.74 3639.12
04)Readmission3 to Readmission4 RP08 11 3604.67 3626.67 3627.08 3675.86
04)Readmission3 to Readmission4 Generalized gamma 4 3619.30 3627.30 3627.37 3645.19
04)Readmission3 to Readmission4 RP09 12 3604.95 3628.95 3629.44 3682.61
04)Readmission3 to Readmission4 Generalized gamma (original) 4 3621.96 3629.96 3630.02 3647.85
04)Readmission3 to Readmission4 RP10 13 3605.14 3631.14 3631.71 3689.28
04)Readmission3 to Readmission4 Gompertz 3 3626.04 3632.04 3632.08 3645.46
04)Readmission3 to Readmission4 Log-logistic 3 3626.47 3632.47 3632.51 3645.89
04)Readmission3 to Readmission4 Exponential 2 3640.60 3644.60 3644.62 3653.54
04)Readmission3 to Readmission4 Weibull (AFT) 3 3639.23 3645.23 3645.27 3658.65
04)Readmission3 to Readmission4 Weibull (PH) 3 3639.23 3645.23 3645.27 3658.65
04)Readmission3 to Readmission4 Gamma 3 3640.18 3646.18 3646.22 3659.60


Session Info

Sys.getenv("R_LIBS_USER")
## [1] "C:/Users/CISS Fondecyt/OneDrive/Documentos/R/win-library/4.0"
rstudioapi::getSourceEditorContext()
## Document Context: 
## - id:        'EA74E2C1'
## - path:      'G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Matching_Process15_JUN_21.Rmd'
## - contents:  <1528 rows>
## Document Selection:
## - [1476, 4] -- [1476, 4]: ''
#save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_carla.RData")

if (grepl("CISS Fondecyt",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/CISS Fondecyt/OneDrive/Escritorio/SUD_CL/mult_state_15_jun.RData")
  } else if (grepl("andre",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("C:/Users/andre/Desktop/SUD_CL/mult_state_15_jun.RData")
  } else if (grepl("E:",rstudioapi::getSourceEditorContext()$path)==T){
    save.image("E:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_jun.RData")
  } else {
    save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/mult_state_15_jun.RData")
  }

sessionInfo()
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Chile.1252  LC_CTYPE=Spanish_Chile.1252   
## [3] LC_MONETARY=Spanish_Chile.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Chile.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gurobi_9.1-0            slam_0.1-47             radiant.update_1.4.1   
##  [4] Rfast_2.0.3             RcppZiggurat_0.1.6      Metrics_0.1.4          
##  [7] muhaz_1.2.6.1           flexsurv_2.0            mstate_0.3.1           
## [10] DiagrammeR_1.0.6.1.9000 Amelia_1.7.6            Rcpp_1.0.5             
## [13] igraph_1.2.6            eha_2.8.1               cobalt_4.2.3           
## [16] MatchIt_3.0.2           tableone_0.12.0         stargazer_5.2.2        
## [19] reshape2_1.4.4          gridExtra_2.3           foreign_0.8-80         
## [22] survMisc_0.5.5          ggfortify_0.4.10        survminer_0.4.8        
## [25] ggpubr_0.4.0            epiR_1.0-15             forcats_0.5.0          
## [28] purrr_0.3.4             readr_1.3.1             tibble_3.0.3           
## [31] tidyverse_1.3.0         dplyr_1.0.1             treemapify_2.5.3       
## [34] sf_0.9-3                ggiraph_0.7.0           finalfit_1.0.1         
## [37] lsmeans_2.30-0          emmeans_1.4.8           RColorBrewer_1.1-2     
## [40] panelr_0.7.3            lme4_1.1-23             Matrix_1.2-18          
## [43] data.table_1.13.0       codebook_0.9.2          devtools_2.3.0         
## [46] usethis_1.6.1           sqldf_0.4-11            RSQLite_2.2.0          
## [49] gsubfn_0.7              proto_1.0.0             broom_0.7.0            
## [52] zoo_1.8-8               rbokeh_0.5.1            janitor_2.0.1          
## [55] plotly_4.9.2.1          kableExtra_1.1.0        Hmisc_4.4-0            
## [58] Formula_1.2-3           survival_3.1-12         lattice_0.20-41        
## [61] ggplot2_3.3.2           stringr_1.4.0           stringi_1.4.6          
## [64] tidyr_1.1.1             knitr_1.29              matrixStats_0.56.0     
## [67] boot_1.3-25            
## 
## loaded via a namespace (and not attached):
##   [1] estimability_1.3    coda_0.19-3         acepack_1.4.1      
##   [4] bit64_0.9-7         multcomp_1.4-13     rpart_4.1-15       
##   [7] generics_0.0.2      callr_3.4.3         TH.data_1.0-10     
##  [10] mice_3.11.0         ggfittext_0.9.0     chron_2.3-55       
##  [13] bit_1.1-15.2        webshot_0.5.2       xml2_1.3.2         
##  [16] lubridate_1.7.9     assertthat_0.2.1    xfun_0.16          
##  [19] hms_0.5.3           evaluate_0.14       fansi_0.4.1        
##  [22] dbplyr_1.4.4        readxl_1.3.1        km.ci_0.5-2        
##  [25] DBI_1.1.0           htmlwidgets_1.5.1   ellipsis_0.3.1     
##  [28] crosstalk_1.1.0.1   backports_1.1.7     insight_0.9.0      
##  [31] survey_4.0          vctrs_0.3.2         remotes_2.2.0      
##  [34] sjlabelled_1.1.6    abind_1.4-5         withr_2.2.0        
##  [37] pryr_0.1.4          checkmate_2.0.0     prettyunits_1.1.1  
##  [40] cluster_2.1.0       lazyeval_0.2.2      crayon_1.3.4       
##  [43] labeling_0.3        pkgconfig_2.0.3     units_0.6-6        
##  [46] nlme_3.1-148        pkgload_1.1.0       nnet_7.3-14        
##  [49] rlang_0.4.7         lifecycle_0.2.0     sandwich_2.5-1     
##  [52] modelr_0.1.8        cellranger_1.1.0    tcltk_4.0.2        
##  [55] rprojroot_1.3-2     KMsurv_0.1-5        carData_3.0-4      
##  [58] reprex_0.3.0        base64enc_0.1-3     processx_3.4.3     
##  [61] png_0.1-7           viridisLite_0.3.0   parameters_0.8.2   
##  [64] KernSmooth_2.23-17  visNetwork_2.0.9    pander_0.6.3       
##  [67] blob_1.2.1          classInt_0.4-3      jpeg_0.1-8.1       
##  [70] rstatix_0.6.0       ggeffects_0.15.1    ggsignif_0.6.0     
##  [73] scales_1.1.1        memoise_1.1.0       magrittr_1.5       
##  [76] plyr_1.8.6          hexbin_1.28.1       compiler_4.0.2     
##  [79] snakecase_0.11.0    cli_2.0.2           ps_1.3.3           
##  [82] htmlTable_2.0.1     MASS_7.3-51.6       tidyselect_1.1.0   
##  [85] highr_0.8           mitools_2.4         jtools_2.0.5       
##  [88] yaml_2.2.1          latticeExtra_0.6-29 grid_4.0.2         
##  [91] tools_4.0.2         parallel_4.0.2      rio_0.5.16         
##  [94] rstudioapi_0.11     uuid_0.1-4          gistr_0.5.0        
##  [97] farver_2.0.3        sjPlot_2.8.4        digest_0.6.25      
## [100] quadprog_1.5-8      car_3.0-8           performance_0.4.8  
## [103] httr_1.4.2          gdtools_0.2.2       effectsize_0.3.2   
## [106] sjstats_0.18.0      colorspace_1.4-1    rvest_0.3.6        
## [109] fs_1.5.0            splines_4.0.2       statmod_1.4.34     
## [112] sessioninfo_1.1.1   systemfonts_0.2.3   xtable_1.8-4       
## [115] jsonlite_1.7.0      nloptr_1.2.2.2      testthat_2.3.2     
## [118] R6_2.4.1            pillar_1.4.6        htmltools_0.5.0    
## [121] glue_1.4.1          minqa_1.2.4         deSolve_1.28       
## [124] class_7.3-17        codetools_0.2-16    maps_3.3.0         
## [127] pkgbuild_1.1.0      mvtnorm_1.1-1       numDeriv_2016.8-1.1
## [130] curl_4.3            BiasedUrn_1.07      zip_2.1.1          
## [133] openxlsx_4.1.5      rmarkdown_2.6       desc_1.2.0         
## [136] munsell_0.5.0       e1071_1.7-3         labelled_2.5.0     
## [139] sjmisc_2.8.5        haven_2.3.1         gtable_0.3.0       
## [142] bayestestR_0.7.2